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Abstract 

In this paper we develop two models for the steady states and evolution of two 
dimensional isothermal self gravitating and rotating incompressible gas which are based 
on the hydrodynamic equations for stratified fluid. The first model is for the steady 
states of the gas while the second addresses the time evolution of the gas subject to 
some constraints. These models reduce the initial five partial differential equations 
that govern this system to two for the steady state model and to three for the time 
dependent model. Analytical and numerical solutions of the model equations are used 
to study the structure of the resulting steady and time dependent states of the fluid 
with some possible astrophysical applications. 


1 Introduction 


The steady states of self gravitating fluid in three dimensions have been studied by a long 
list of illustrious mathematical physicists. (For an extensive list of references see [1,2,3]). 
The motivation for this research was clue to the interest in the shape, stability and evolution 
of celestial bodies and systems[12], We now know however that many celestial objects such 
as galaxies and our solar system exhibit (effectively) ’’two dimensional structure” [4,5,6,7]. 
Furthermore recent discoveries are leading us to believe that systems similar to our solar 
system are ” abundant” in the galaxy and their existence might be due to the collapse of 
a two dimensional interstellar cloud under gravitation (this is the so called the ’’nebular 
theory”) [8,9,10,11,17]. This data leads us to believe that there is a fundamental physical 
process which we do not understand fully as yet that leads to the formation of planetary 
systems throughout the galaxy (and beyond). 

This background motivates us to investigate in this paper the steady states and time- 
dependent evolution of a self gravitating and rotating fluid in two dimensions. This problem 
has been explored by a large number of investigators using elaborate analytic methods and 
computer simulations which involve, in general, thermodynamic considerations, magnetohy¬ 
drodynamics modeling and turbulence.(For a complete list of references see [8,9,10,13,16]). 
While these are important issues we still need, in our opinion, prototype analytic models 
that are able to capture the evolution of this process and lead to insights about its possible 
outcomes. 

In this paper we attempt to develop such a model using the basic hydrodynamic equations 
that govern the time-dependent evolution of an isothermal, incompressible, stratified (i.e 
non constant density) and rotating fluid in two dimensions under gravity [1,2,3]. (The 
justification for the reduction from three to two dimensions has been discussed by many 
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authors. A lucid treatment is given in Ref. [12] pp.1-12).Under these assumptions we show 
that the number of model equations can be reduced from five to a system of two equations 
for the steady states and three coupled equations for its evolution. The models contain 
some ” parameter functions” which encode information about the asymptotic mass density 
distribution of the fluid and its momentum. 

The steady state model was investigated by us previously. [18,19,20] However in this paper 
we consider a more general model in which the ”gas cloud” is rotating also with uniform 
angular velocity c o and study the impact of this rotation on the matter distribution in the 
steady state. 

To study the predictions of these models we use both analytical and numerical methods to 
solve their equations under a variety of conditions. In particular we consider radial solutions 
to these equations which represent the evolution of an interstellar cloud with isothermal 
equation of state [10]. 

It might be argued that the hydrodynamic assumptions we are making in this paper are 
not realistic from astrophysical point of view. However our main goal is to capture analyt¬ 
ically, as far as possible, the nonlinear and time dependent aspects of the processes under 
consideration. Accordingly our results might be useful to provide some analytic insights and 
guidelines for more elaborate work on this topic. 

The plan of the paper is as follows: In Sec 2 we present the basic hydrodynamic equations 
and show how one can reduce them to a coupled system of three equations. Sec 3 presents 
further simplifications of these equations. The first is for the steady states of the model. The 
second is for the time dependent evolution of the gas cloud under the assumption of constant 
vorticity. In Sec 4 we present analytical and numerical radial solutions of these equations. 
We end up Sec 5 with summary and conclusions. 
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2 Derivation of the Model Equations 


Following the standard convention [1,14,15] we model the time dependent non-relativistic 
flow of an incompressible fluid in two dimensions (x,y) by the hydrodynamic equations of 


inviscid and incompressible stratified fluid 

u x + Vy = 0 (2.1) 

Pt T Ujp x V Py 0 (2.2) 

put + p{uu x + vUy) = -p x - P(p x + poj 2 x (2.3) 

pU t + p(UV X + Wy) = ~Py - p(j)y + pij ]) (2.4) 

v 2 0 = AitGp (2.5) 


where subscripts indicate differentiation with respect to the indicated variable, u = {u,v) 
is the fluid velocity, p is its density, p is the pressure, (j) is the gravitational field and G is 
the gravitational constant. The terms pu 2 x , pu 2 y represent the components of the apparent 
centrifugal force due to the rotation of the gas cloud with angular velocity co. 

We can nondimensionalize these equations by introducing the following scalings 

t = jj-, x = Lx, y = Ly, u = U 0 u , v = U 0 v, p = p 0 p, p = poU^p, 0 = £/ o 2 0, u = 

Uq L 

( 2 . 6 ) 

where L, Uq, po are some characteristic length,velocity and mass density respectively that 
characterize the problem at hand. Substituting these scalings in eqs. (12.If) - (12.51) and drop¬ 
ping the tildes these equations remain unchanged (but the quantities that appear in these 
equations become nondimensional) while G is replaced by G = Gp ^ . (Once again we drop 
the tilde). 

In view of eq. (12.ip we can introduce a stream function so that 

u = Vv, U = -Vtr • (2.7) 
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Using this stream function we can rewrite eq. (|2.2[) as [13,15] 


Pt + J{p,^} = 0 

(2.8) 

where for any two (smooth) functions /, g 


, , df dg df dg 

U ’^ 1 Ur &y dy dx 

(2.9) 

Using ib the momentum eauations (12.31).(12.41) become 


PipPyt + Ipy^yx ~ ^x^yy) = ~Px ~ pfx + PW 2 £ 

(2.10) 

P(-Vht - V’yV’x* + ^xf>xy) = ~Py ~ Pfy + P^V 

(2.11) 


To eliminate p from these equations we differentiate eq. (I2.10p and eq. (12. lip with respect 
to y,x respectively and subtract. This leads to 

Pyi^Pyt T Ipylpyx Ipx'lpyy) T Pfyyyt T 4’ y^Pyyx ^x^Pyyy) (2.12) 

Px ( Wxt 'Ipy'lpxx T ’Ipx'lpxy) P( *Pxxt '^Py'^Pxxx T 'Ipx'^Pxxy) P\ T f i P\ 

where r 2 = x 2 +y 2 . The sum of the second and fourth terms in this equation can be rewritten 
as 

p(V 2 ^ + pJ{V 2 ^}- (2.13) 

To reduce the first and third terms in (12.121) we use (12.8p . We obtain 

P!J ( 4'yt T PppPyx Uxht/y) Px ( Ux/. PppPxx T UxUxy) (2.14) 

PyipPyt 4~ fp/PxyPpjx (pt T PxPpy)Ppjy T PxUx/ 4“ (PpxPy Pt ) Uxx PxpPxPpxiy 

Py^yt + Py^yt ~~ Pt^ 2 ip + ^</{(Vh) 2 + (^) 2 , P>- 
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combining the results of (j2. 13[) and (12.14)1 eq. (12.121) becomes 


py'ijjyt + p x ip xt - + p(V 2 00 + pJ{V 2 0, 0} + i J{(0x) 2 + (0j/) 2 , p} (2.15) 

= + J{\^ 2 ,p) 

Thus we have reduced the original five equations (I2.ip - (I2.5I) to three equations (12.5p . (12.8p 
and (12.151) . Although (I2.15P is rather cumbersome in general, it can simplified further under 
some restrictions which are presented in the following section. 

3 Simplification of the Model Equations 

Equation (I2.15P can be simplified further in two cases. The first is when we consider only 
steady states of the flow and the second is when the flow vorticity is constant. 

3.1 A Model for the Steady States 

When we consider only steady states of the flow (12 ,8j) implies that 0 = 0(p) and after some 
algebra [18] (12.15p reduces to 

H(p)^S7iH( P y/ 2 S7p) +t- Av = S(p). (3,1) 

Where 

H(p) = ptPl (3.2) 

and S(p) is some function of p. Thus the equations governing the steady state are (13.ip . 
(12.5ft and H(p) and S(p) are ’’parameter functions” which determine the nature of the steady 
state. 
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3.1.1 The Physical Meaning of the Functions H(p), S(p ) 


The function H(p) is a parameter function which is determined by the momentum (and 
angular momentum) distribution in the fluid. From a practical point of view the choice of this 
function determines the structure of the steady state density distribution. The corresponding 
flow field can be computed then aposteriori (that is after solving for p) from the following 
relations. [18] 


u = 


H{p)dp 


v = 


'H(p) dp 
p dx 


(3.3) 


p dy ’ 

The function S(p) that appears in eq. (13.11) can be determined from the asymptotic 
values of p and (j) on the boundaries of the domain on which eqs (I2.5p . fl3.ip are solved. When 
these asymptotic values are imposed or known one can evaluate the left hand side of eq. 
(13. 1 1) on the domain boundaries and re-express it in terms of p only to determine S(p) on 
the boundary of the domain. However, the resulting functional relationship of S on p must 
then hold also within the domain itself since S does not depend on x, y directly. 

For example if we assume that on an infinite domain h(p) — 1, u — 0 and the asymptotic 
behavior of p and (j) is given by 


lim p(r) = e Qr2 , lim <f>(r) = 4 a 2 r 2 e ar2 (3.4) 

r—»• oo r—too 

then (asymptotically) (13. ip evaluates to 

S(p) = —4 ae~ ar = — Aap (3.5) 


3.2 A model for the Time Evolution 

To begin with we consider the case where the vorticity is zero and then generalize to the 
case where the flow vorticity is constant. 
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When the flow vorticity Vxuis zero then V 2, 0 = 0 eq. (12.15j> becomes 


Py^yt + Px'tpxt + y) 2 ,P} = -J{(/>,p} (3.6) 

However when the vorticity is zero we can introduce the velocity potential r] which satisfies 
rj x = u, r]y = v. Replacing ^ by r/ in (13. 6 p we obtain 

J{Vt + ^ [(Vx) 2 + {Vy) 2 } + 0) p} = 0 (3.7) 

Hence 

Vt + ^\(Vx) 2 + {Vy) 2 ] + (p = S(p) (3.8) 

The equations of the flow in this case are 

Pt + VxPx + VyPy = o (3.9) 

(which replaces (12. 8 p ). (13.8j) and (12.5ft . 

To generalize this reduction to the case where V 2, 0 = a (where a is any constant) we 
define 

vi = ip y , v 2 = -iIj x + ax 

Therefore 

{vi) y - ( v 2 ) x = 0 , 

which implies that there exists a function 77 so that 


Hence 


Vx = v 1 , rjy — v 2 . 


r) x = ^ y , Vy = -'fpx + ax 


( 3 . 10 ) 










Using these relations to substitute rj for 0 in (13.121) leads to 


PyVxt Px (jly ax)t T 


-ap t 




ax) 2 + ( p x ) 2 ,p } 


= -J{0,p}- 


(3.11) 


Therefore 


J{Vt,p} ~ ap t + -J{(v y 


ax) 2 + {pxfip} 


-J{(p,p}. 


(3.12) 


Hence 


- ap t + J{p t + -[(% 


ax) 2 + (7i x ) 2 } + 0,p} 


0 . 


(3.13) 


Using (12. 8 p we have 


"■/{t-'./'l + J{ri, + -[(i), - ax) 2 + (ife) 2 ] + cp,fi} 


(3.14) 


It follows then that 


- 00 + Vt+ xliVv ~ ax) 2 + (rfe) 2 ] + <p = S(p). 


(3.15) 


If a Y 0, 0 can be eliminated from this equation if we differentiate with respect to y and use 
(13.1 Op to obtain 


- arj x + 


Vt + 7 ,[(Vy - ax) 2 + (px) 2 } + 0 


= s( P \ 


(3.16) 


4 Radial Solutions for the Steady State Model 

When we consider the special case where in polar coordinates p = p(r) and 0 = 0(r) the 
system (13.ip and (12. 5 p with H(p) = 1 reduces to 

/ = -^ + S(p)-*+;hV (4.1) 

r 2 

0" = —0' + cp, c = 47 tG (4.2) 

r 
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To solve this system of equations we let S(p) = ap , solve (14.ip for 0 and substitute the result 
in (14.211 . This leads to the following fourth order equation for p 

p"" H— p>" — fa + —\ p" + f— \ p' + cp = 2u 2 . (4.3) 

rtn l ir>i J l y»0 rtr> J 

The general solution of this equation is 
2u 2 

P = -b C'iJo(air) + C 2 J 0 (bir) + C 3 Y 0 (a 1 r) + C A Y Q (bir) (4.4) 

c 

where Jo and Y 0 are Bessel functions of the first and second kind of order 0 and 

ci\ = - y—2a + 2 \Ja 2 — 4c + ofi, b\ = -\j—2a — 2y/a 2 — 4c + 

Assuming no singularity at the origin we set C 3 = C 4 = 0. To assess the impact of the 
rotation term on the steady state we solved this system for C \, C 2 on a circular disk using 
the boundary conditions p(0) = 1 and p( 8 ) = 0 with c = 1, a = —19.4. The results of these 
computations for different values of u are plotted in Fig. 1. In this figure we see that the 
separation between the density peaks become more pronounced as 00 increases. This might 
interpreted as leading to the creation of protoplanets around the central core. 

A strong dependence on cu is shown in Fig. 2 which has the same parameters as Fig. 
1 except that the boundary conditions on p are: p(0) = 0.35 and p( 8 ) = 0.25. This figure 
illustrate clearly the effect that rotation can have on the pattern of density fluctuations within 
the cloud. Furthermore in this figure the magnitude of the density fluctuations reverses itself 
as 00 becomes larger viz. the higher density peaks are placed at larger values of r. (Which is 
reminiscent of the situation in the solar system) 
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5 Radial Solutions for the Time Evolution Model 


The system (I2.5j) . fl3.8j) and (13.91) can be simplified further if we use polar coordinates and 
assume that p, r), (ft are functions of r and t only. We obtain, 


(5.5) 


Pt T PrPr 0, 


(ftrr T — ( ftr — Cp — 0? 


r 



where c = 47 tG. 

5.1 Steady States 

When we consider a steady state solutions of (15.51) then p t = 0 and p t = 0. If follows from 
the first equation in (j5.5[) that either p r or must be zero. In the first case p is constant 
and we can let p = 1 without loss of generality. When r/ r is zero we must have (ft = S(p) and 
the second equation in (15.51) becomes 



(5.6) 


where primes denote differentiation with respect to p. 

We consider these two cases separately. 

A. Steady state with p = 1 

Since p = 1 the function S(p) is a constant and the general solution for (ft is 


(ft = C -r 2 + Ci In r + C 2 - 


(5.7) 


where Cf, C 2 are arbitrary constants. The equation for p becomes 


= S- -C - Ci ln(r) - C 2 . 


(5.8) 
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(S can be absorbed in C 2 but we leave it in this form as these two constants have different 
physical meaning). If we let C\ — 0 to avoid the singularity at the origin (j5.8[) yields 


rj = ± 


-r\/8S~ 


8 C 9 — 2 cr 2 + —— arc tan 

\/c 


\l2ci 


V 8 S - 8 C 2 - 2cr 2 


+ C 3 (5.9) 


B. Steady states with rj — 1 

In this case the solution of (15. 6 p depends on the nature of the function S(p). In general 
this equation has to be solved numerically. However we present here analytical solutions of 
this equation for two special cases. 


1. S(p ) = ap where a is a constant. The solution to (15.6p in this case is 

p = C 4 J 0 (y^r) +C 5 V 0 (y^r) (5.10) 

It follows then that the nature of the steady state is determined by the ratio A 
sample of the resulting p profiles is presented in Fig. 3. To obtain this figure we 
considered a pinched disk with p(0.01) = 1, p'(0.01) = —10 and c = a. The resulting 
steady state has an increase in the material density towards the circumference of the 
disk. Similar graphs were obtained numerically for S(p) = ap n , n = 2 , 3. 


2 . S(p) = a ln(p) 


In this case we have 


where 


P = 


2C\cr 2 cos( 6>) 2 


(5.11) 


6 = 


- (In r - C 2 ) 


2^a 

Substituting Ci = C 2 = c = a = 1 we obtain Fig. 4 which might be interpreted as 
representing a binary system. 
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5.2 Perturbations from the steady state p = 1 


We consider in this section a disk of radius 1 with a steady state po = 1 and S(p) = 0. 
Letting 0(1) = 0 and using (15.71) . (15.91) (with Cj = 0) this yields the following equations for 
the steady state 


Vo 


M r ) = ^(r 2 -1) 


y/2c 


arcsin(r) + r\f[ 1 — r 2 ) 


For a perturbation from this state, viz. 


p(t,r) = po + epi(«,r), 0(i, r) = <#> 0 + r), i?(i, r) = jfa + ei)i(i,r) 


(5.12) 


we obtain to first order in e the following system of equations: 


(Pi)f H - 2-(P 1 )?- ~ 

(0l)rr “1“ — (0l)r 0, 




vW 


-(//l) r + 0i = o. 


(5.13) 


The equation for p 1 in (j5. 13[) can be solved analytically. Its general solution is 


pi = F ( \j - arcsin(r) — t 


(5.14) 


where F is any smooth function of its variable which has to be adjusted to the initial 
conditions of the perturbation. The second equation in (15.131) is a (reduced) Poisson equation 
and its general solution can be expressed by quadratures 

arcsin(r) — tj dr + F x (t)} 


>1 = 


dr + F 2 (t ) 


(5.15) 


where F x (t), F 2 (t) have to be determined by the boundary conditions on 0]. Finally one can 
obtain also an expression for the solution for r/ in terms of quadratures. 
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For example if the initial perturbation in p is pi(0,r) = ar where a is a constant then 


and 



Pi(t,r) = a sin ^arcsin(r) — 


= a < r cos 


f]-vT= 


r z sm 


(5.16) 


The evaluation of 0! using the second equation in (I5.13p and (15.1611 is straightforward. It 
should be obvious how one can generalize this example to other expressions for pi(0 ,r). 

A second approach to the solution of the system (15.131) is to assume exponential depen¬ 
dence in time, viz. 

pi = e at R(r ), rji = e at E(r), 0i = e at P(r) (5.17) 


This ansatz reduces (15.13jl to a system of ordinary differential equations 

V / 2 C (^ ~ r2 ) R{r)' + aR(r) = 0, (5.18) 

P(r)" + rP{r)' — crR{r ) = 0, 
v/2c(l ~ r2 ) E( r y + a E(r) + P(r) = 0. 


As before the equation for R(r ) can be solved analytically, 


R(r ) = C\ exp ( —a\j - arcsin(r) 

while the equations for P(r) and E(r) can be solved by quadratures or numerically. 
A numerical approach to the solution for 0! and r)\ in (15.13[) is also possible. 


5.3 Perturbations from the steady state r )o = 1 

As in the previous subsection we consider again a disk of radius 1 and let S(p ) = ap. The 
general steady state solution for p is given by (I5.9I) . Assuming no singularities in p (ie. no 
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protostar at the origin) we must set C 2 = 0 in this equation. Furthermore since p > 0 it 
follows that we must have yj—c/a = f3 where f3 is the first zero of Jq. Thus 


Po = MM, 0o = «Po 

(where we normalized rho 0 at r = 0 to be 1). 

For a perturbation from this steady state in the form given by (I5.12p we obtain to first 
order in e the following system of equations: 

(0i)tt + — (0i)r - cpi = 0, (5.19) 

r 

(rji )t ~ & pi + 0i = 0 

(Pl)t - PMM(Vl)r = 0 

Where J\ is Bessel function of the first kind of order 1. The evolution of an initial pertur¬ 
bation pi = exp(—5r ) with a = 0.01 from the steady state is plotted in Fig5. This figure 
shows that as time progresses there is an accumulation of matter near the center of the disk. 
At the same time there is an initial separation between the core and the rest of the disk. 

We computed also the solution to the system (15.51) with an initial matter distribution 
p(0,r) = 1+S11 k 27rr ) ; c — o.Ol and t G [0,9]. The results of the simulation (Fig. 6) show that 
as time progresses matter is starting to build up in the vicinity of the center of the disk and 
around r — 1. At the same time there is a decrease in matter density in between these two 
points. 

6 Summary and Conclusions 

In previous publications [18-20] we treated only the steady states of two dimensional self 
gravitating fluid. In this paper we generalized this model to include disk rotation and 
assessed the impact of this addition on the distribution of matter in the disk. We were able 
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also to address the time dependent evolution of this fluid under restrictions on its vorticity. 
This enabled us to simplify considerably the equations which govern its evolution. While this 
is a highly idealized model in the context of astrophysical applications it may still provide 
some analytical insights for more elaborate models. 

In this paper we considered only radial solutions of this model. More general solutions 
which are not radial will have to be explored next. 
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List of Captions 


Fig. 1 Steady states with a = —19.4, c = 1 and boudary conditions p(0) = 1, 
p(8) = 0 with different values of u> 

Fig. 2 Steady states with a = —19.4, c = 1 and boudary conditions p(0) = 0.35, 

p(8) = 0.25 with different values of oj 

Fig. 3 The steady state that corresponds to (15. 10(1 

Fig. 4 The steady state that corresponds to (15.111) 

Fig. 5 Using (15. 19ft to solve for pi with a = —0.01, 

c = 0.0578 and initial perturbation pi = exp(—5r) with a protostar at the origin. 
Fig. 6 Using (15.51) to compute the evolution of p. The initial matter distribution is 
p(0,r) = 1 +sl U 27rr ) anc [ c — o.Ol. No protostar at the origin. 
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Evolution of initial perurbation in p=exp(-5*r) 



Figure 5: 












Time Evolution of p with p(0,r)=(1 +sin2jt r))/2 



Figure 6: 












